# Denmark - Municipalities

#libraries
library(tidyverse) 
library(Hmisc)
library(lubridate)
library(broom)
library(gganimate)
library(transformr)
library(gtools)
library(readxl)
library(stringr)
library(estimatr)
library(texreg)

# loading the google data
dk_local <- read.csv("Global_Mobility_Report_August11.csv")

dk_local2 <- dk_local %>% 
  filter(country_region=="Denmark")

dk_local2$sub_region_2 <- dk_local2$sub_region_2 %>% 
  str_replace(" Municipality", "")

dk_local2 <- dk_local2 %>% 
  mutate(sub_region_2 = case_when(sub_region_2=="Copenhagen" ~ "København",
                                  sub_region_2=="Brondby" ~ "Brøndby",
                                  sub_region_2=="Nordfyn" ~ "Nordfyns",
                                  sub_region_2=="Vesthimmerland" ~ "Vesthimmerlands",
                                  TRUE ~ sub_region_2))

# loading vote shares for Danish People's Party
DPP <- read_excel("Denmark background data/danish_peoples_party.xlsx",
                  range="FVPANDEL!B5:C103",col_names = FALSE, na = ":")
names(DPP) <- c("sub_region_2","vote_share_dpp")

# loading vote shares for Nye Borgerlige
NB <- read_excel("Denmark background data/nye_borgerlige.xlsx",
                 range="FVPANDEL!B5:C103",col_names = FALSE, na = ":")
names(NB) <- c("sub_region_2","vote_share_nb")

# loading vote shares for Stram Kurs
SK <- read_excel("Denmark background data/stram_kurs.xlsx",
                 range="FVPANDEL!B5:C103",col_names = FALSE, na = ":")
names(SK) <- c("sub_region_2","vote_share_sk")

# loading vote shares for Social Democrats
SD <- read_excel("Denmark background data/social_democrats.xlsx",
                 range="FVPANDEL!B5:C103",col_names = FALSE, na = ":")
names(SD) <- c("sub_region_2","vote_share_sd")

# house prices (fourth quarter 2019)
house_prices <- read_excel("Denmark background data/house_prices.xlsx",
                           range="BM010!C4:D101",col_names = FALSE, na = ":")
names(house_prices) <- c("sub_region_2","house_prices")


# Density
density <- read_csv("Denmark background data/density.csv")

density$sub_region_2 <- density$sub_region_2 %>% 
  str_replace(" Kommune", "")

density <- density %>% 
  mutate(sub_region_2 = case_when(sub_region_2=="Københavns" ~ "København",
                                  sub_region_2=="Bornholms" ~ "Bornholm",
                                  TRUE ~ sub_region_2))

# post secondary education
post_secondary <- read_csv("Denmark background data/post_secondary.csv")

post_secondary$sub_region_2 <- post_secondary$sub_region_2 %>% 
  str_replace(" Kommune", "")

post_secondary <- post_secondary %>% 
  mutate(sub_region_2 = case_when(sub_region_2=="Københavns" ~ "København",
                                  sub_region_2=="Bornholms" ~ "Bornholm",
                                  TRUE ~ sub_region_2))

# lower secondary education
lower_secondary <- read_csv("Denmark background data/lower_secondary.csv")

lower_secondary$sub_region_2 <- lower_secondary$sub_region_2 %>% 
  str_replace(" Kommune", "")

lower_secondary <- lower_secondary %>% 
  mutate(sub_region_2 = case_when(sub_region_2=="Københavns" ~ "København",
                                  sub_region_2=="Bornholms" ~ "Bornholm",
                                  TRUE ~ sub_region_2))

# full-time unemployed
unemp <- read_csv("Denmark background data/full_time_unemployed.csv")

unemp$sub_region_2 <- unemp$sub_region_2 %>% 
  str_replace(" Kommune", "")

unemp <- unemp %>% 
  mutate(sub_region_2 = case_when(sub_region_2=="Københavns" ~ "København",
                                  sub_region_2=="Bornholms" ~ "Bornholm",
                                  TRUE ~ sub_region_2))


# Average disposable income
income <- read_excel("Denmark background data/disposable_income.xlsx",
                     range="INDKP106!E4:F101",col_names = FALSE, na = ":")
names(income) <- c("sub_region_2","disp_inc")

# Net wealth
net_wealth <- read_excel("Denmark background data/net_wealth.xlsx",
                         range="FORMUE2!B4:C101",col_names = FALSE, na = ":")
names(net_wealth) <- c("sub_region_2","net_wealth")

# Pop over 65 data
pop65 <- read_csv("Denmark background data/plus65.csv")

pop65$sub_region_2 <- pop65$sub_region_2 %>% 
  str_replace(" Kommune", "")

pop65 <- pop65 %>% 
  mutate(sub_region_2 = case_when(sub_region_2=="Københavns" ~ "København",
                                  sub_region_2=="Bornholms" ~ "Bornholm",
                                  TRUE ~ sub_region_2))

# Pop under 18 data
pop18 <- read_csv("Denmark background data/under18.csv")

pop18$sub_region_2 <- pop18$sub_region_2 %>% 
  str_replace(" Kommune", "")

pop18 <- pop18 %>% 
  mutate(sub_region_2 = case_when(sub_region_2=="Københavns" ~ "København",
                                  sub_region_2=="Bornholms" ~ "Bornholm",
                                  TRUE ~ sub_region_2))


# Employment by sector
sector <- read_excel("Denmark background data/employment_by_sector.xlsx",
                     range="RAS301!A5:AL104",col_names = FALSE, na = ":")
names(sector) <- c("sub_region_2","A","B","CA","CB","CC","CD","CE","CF","CG","CH","CI","CJ","CK","CL","CM", 
                   "D", "E", "F", "G", "H", "I", "JA", "JB", "JC", "K", "L", "MA", "MB", "MC", "N", "O", "P", "QA", "QB", "R", "S", "X")
sector$C <- sector$CA+sector$CB+sector$CC+sector$CD+sector$CE+sector$CF+sector$CG+sector$CH+sector$CI+sector$CJ+sector$CK+sector$CL+sector$CM
sector$J <- sector$JA+sector$JB+sector$JC
sector$M <- sector$MA+sector$MB+sector$MC
sector$Q <- sector$QA+sector$QB
sector$total_employment <- sector$A+sector$B+sector$C+sector$D+sector$E+sector$F+
  sector$G+sector$H+sector$I+sector$J+sector$K+sector$L+sector$M+sector$N+sector$O+
  sector$P+sector$Q+sector$R+sector$S+sector$X

# farming, forestry, fishing
sector$emp_farm_fish_forest <- sector$A/sector$total_employment

# industry
sector$emp_industry <- (sector$B+sector$C+sector$D+sector$E)/sector$total_employment
sector$emp_construction <- sector$F/sector$total_employment

# wholesale, retail, transport, hotel, restaurants
sector$emp_wholesale_retail_transport_acc <- (sector$G+sector$H+sector$I)/sector$total_employment

# financial services, real estate, profesionals, scientific, admin
sector$emp_finance_lib <- (sector$K+sector$L+sector$M+sector$N)/sector$total_employment

# frontline workers: defense, police, fire, educational sector, healthcare sector
sector$emp_adm_health_education <- (sector$O+sector$P+sector$Q)/sector$total_employment

pr_sector <- subset(sector, select = c("sub_region_2",
                                       "emp_finance_lib","emp_construction","emp_wholesale_retail_transport_acc","emp_adm_health_education"))

# merging the data frames
DK<-pop65 %>% 
  left_join(DPP, by  = "sub_region_2") %>% 
  left_join(NB, by  = "sub_region_2") %>% 
  left_join(SK, by  = "sub_region_2") %>% 
  left_join(SD, by  = "sub_region_2") %>% 
  left_join(house_prices, by  = "sub_region_2") %>% 
  left_join(net_wealth, by  = "sub_region_2") %>% 
  left_join(income, by  = "sub_region_2") %>% 
  left_join(post_secondary, by  = "sub_region_2") %>% 
  left_join(lower_secondary, by  = "sub_region_2") %>% 
  left_join(density, by  = "sub_region_2") %>% 
  left_join(unemp, by  = "sub_region_2") %>% 
  left_join(pop18, by  = "sub_region_2")  %>% 
  left_join(pr_sector, by  = "sub_region_2") 

DK$vote_share_populists <- DK$vote_share_dpp+DK$vote_share_nb+DK$vote_share_sk


# merging with google data
dk_local2<-dk_local2 %>% 
  left_join(DK, by  = "sub_region_2")

dk_local2 <- dk_local2 %>% 
  filter(sub_region_2!="")

dk_local2  <- dk_local2 %>% 
  mutate(date = as.Date(date)) %>% 
  rename (grocery = grocery_and_pharmacy_percent_change_from_baseline,
          parks = parks_percent_change_from_baseline,
          residential = residential_percent_change_from_baseline,
          retail_rec = retail_and_recreation_percent_change_from_baseline,
          transit = transit_stations_percent_change_from_baseline,
          workplace = workplaces_percent_change_from_baseline,
          location = sub_region_1)

weekend_df <- data.frame(date = unique(dk_local2$date)) %>%
  mutate(day = weekdays(date)) %>%
  mutate(weekend= case_when(
    day %in% c("Saturday", "Sunday") ~ 1,
    date==as.Date("2020-04-09") ~ 1,
    date==as.Date("2020-04-10") ~ 1,
    date==as.Date("2020-04-13") ~ 1,
    date==as.Date("2020-05-01") ~ 1,
    date==as.Date("2020-05-08") ~ 1,
    date==as.Date("2020-05-21") ~ 1,
    date==as.Date("2020-05-22") ~ 1,
    date==as.Date("2020-06-01") ~ 1,
    date==as.Date("2020-06-05") ~ 1,
    date==as.Date("2020-08-22") ~ 1,
    date==as.Date("2020-08-23") ~ 1,
    date==as.Date("2020-08-29") ~ 1,
    date==as.Date("2020-08-30") ~ 1,
    date==as.Date("2020-09-05") ~ 1,
    date==as.Date("2020-09-06") ~ 1,
    TRUE ~ 0))

dk_local2$log_house_prices <- log(dk_local2$house_prices)
dk_local2$log_density <- log(dk_local2$density)
dk_local2$log_net_wealth <- log(dk_local2$net_wealth)


# number of tests pr municipality pr day
dk_tests <- read.csv("Denmark background data/Data-Epidemiologiske-Rapport-01122020-2kma/Municipality_tested_persons_time_series.csv", sep=";")
# from wide to long
dk_tests_long <- gather(dk_tests, sub_region_2, tests, Copenhagen:NA., factor_key=TRUE)
names(dk_tests_long)[names(dk_tests_long)=="PrDate_adjusted"] <- "date"

# number of positive tests pr municipality pr day
dk_confirmed <- read.csv("Denmark background data/Data-Epidemiologiske-Rapport-01122020-2kma/Municipality_cases_time_series.csv", sep=";")
# from wide to long
dk_confirmed_long <- gather(dk_confirmed, sub_region_2, pos_test, Roskilde:Samsø, factor_key=TRUE)
names(dk_confirmed_long)[names(dk_confirmed_long)=="date_sample"] <- "date"

dk_pos_test<-dk_tests_long %>% 
  left_join(dk_confirmed_long, by  = c("sub_region_2", "date"))

dk_pos_test$pos_test_ratio <- dk_pos_test$pos_test/dk_pos_test$tests
dk_pos_test$sub_region_2 <- as.character(dk_pos_test$sub_region_2)

dk_pos_test <- dk_pos_test %>% 
  mutate(sub_region_2 = case_when(sub_region_2=="Copenhagen" ~ "København",
                                  sub_region_2=="Høje.Taastrup" ~ "Høje-Taastrup",
                                  sub_region_2=="Ringkøbing.Skjern" ~ "Ringkøbing-Skjern",
                                  sub_region_2=="Faaborg.Midtfyn" ~ "Faaborg-Midtfyn",
                                  sub_region_2=="Ikast.Brande" ~ "Ikast-Brande",
                                  sub_region_2=="Lyngby.Taarbæk" ~ "Lyngby-Taarbæk",
                                  TRUE ~ sub_region_2))

dk_pos_test$date <- as.Date(dk_pos_test$date)
dk_pos_test$pos_test_ratio[is.nan(dk_pos_test$pos_test_ratio)] <- NA

dk_local2<-dk_local2 %>% 
  left_join(dk_pos_test, by  = c("sub_region_2", "date"))


## Figure 5
# The Danish Dragon
median(dk_local2$vote_share_dpp, na.rm=T)
median(dk_local2$vote_share_sd, na.rm=T)
median(dk_local2$vote_share_nb, na.rm=T)
median(dk_local2$vote_share_sk, na.rm=T)

# DPP
dk_local2 %>% 
  #left_join(weekend_df) %>%
  #filter(weekend==0)  %>%
  ggplot(aes(x = date, y = workplace, group=sub_region_2, alpha = 0.2))+
  geom_jitter(aes(color=as.factor(vote_share_dpp>10.1)))+
  scale_color_manual(values=c('gray','red'))+
  geom_vline(xintercept =as.Date("2020-03-11"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-04-14"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-05-11"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-06-08"), linetype="dashed", color = "blue")+
  ylab("Change in Workplace Activity from Baseline")+xlab("Date")+
  theme_classic()+
  annotate(geom="text", x=as.Date("2020-03-04"), y=60, label="Lockdown", color="red", size = 2.5)+
  annotate(geom="text", x=as.Date("2020-04-01"), y=60, label="Phase 1 of reopening", color="blue", size = 2.5)+
  annotate(geom="text", x=as.Date("2020-04-28"), y=60, label="Phase 2 of reopening", color="blue", size = 2.5)+
  annotate(geom="text", x=as.Date("2020-05-26"), y=60, label="Phase 3 of reopening", color="blue", size = 2.5)+
  theme(legend.position = "none")





## Figure 6
## coef plot by day
# workplace - DPP+NB+SK
coefs_map_workplace <- map_df(unique(dk_local2$date), function(date_i) {
  filtered_data <- filter(dk_local2, date == date_i) %>%
    lm(data =., workplace ~ vote_share_populists+log_house_prices+disp_inc+log_density+emp_finance_lib+emp_construction+emp_wholesale_retail_transport_acc+emp_adm_health_education+plus65+under18+I(location)) %>%
    broom::tidy() %>%
    rename(coef = estimate, se = std.error) %>%
    mutate(date = date_i)
})

coefs_map_workplace <- coefs_map_workplace %>% 
  mutate(upper = coef+2*se,
         lower = coef-2*se) %>% 
  filter(term %in% c("vote_share_populists","log_house_prices", "disp_inc", "log_density", "under18", "plus65")) %>% 
  mutate(term = case_when(term=="vote_share_populists" ~ "Populist Vote",
                          term=="log_house_prices" ~ "Log House Prices",
                          term=="disp_inc" ~ "Disposable Income",
                          term=="log_density" ~ "Log Density",
                          term == "under18" ~ "Population < 18",
                          term == "plus65" ~ "Population > 65"))

coefs_map_workplace %>% 
  left_join(weekend_df) %>%
  filter(weekend==0)  %>%
  ggplot(aes(x = date, y = coef))+
  geom_line()+
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
  geom_vline(xintercept =as.Date("2020-03-11"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-04-14"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-05-11"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-06-08"), linetype="dashed", color = "blue")+
  geom_ribbon(aes(ymin=lower, ymax = upper, alpha=0.3), color="lightgray")+
  #lockdown
  xlab("Date")+ylab("Day by Day Coefficients (DV: Workplace Activity)")+
  theme_minimal()+
  theme(legend.position = "none") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_wrap(~term, scales = "free_y")




## pooled regressions
# date dummies and lags
dk_local2 <-dk_local2  %>% 
  mutate(after_march_10 = as.numeric(date>as.Date("2020-03-10")),
         after_april_13 = as.numeric(date>as.Date("2020-04-13")),
         after_may_10 = as.numeric(date>as.Date("2020-05-10")),
         after_june_7 = as.numeric(date>as.Date("2020-06-07"))) %>% 
  group_by(sub_region_2) %>% 
  mutate(lag_workplace = lag(workplace),
         l2_workplace = lag(workplace, 2),
         l3_workplace = lag(workplace, 3),
         l4_workplace = lag(workplace, 4),
         l5_workplace = lag(workplace, 5),
         l6_workplace = lag(workplace, 6),
         l6_workplace = lag(workplace, 7)) %>% 
  ungroup()


# removing weekends and national holidays
dk_local2<-dk_local2 %>% 
  left_join(weekend_df, by  = "date")
dk_local2_no_weekend <- dk_local2 %>% 
  filter(weekend==0)



## Table 3
m1 <- lm_robust(data =dk_local2_no_weekend, workplace ~ 
                  lag_workplace+
                  vote_share_populists+after_march_10+after_april_13+after_may_10+after_june_7+
                  after_march_10:vote_share_populists+
                  after_april_13:vote_share_populists+
                  after_may_10:vote_share_populists+
                  after_june_7:vote_share_populists+
                  log_house_prices+disp_inc+log_density+emp_finance_lib+emp_construction+emp_wholesale_retail_transport_acc+emp_adm_health_education+plus65+under18,
                fixed_effects = location, se_type = "HC1")

m2 <- lm_robust(data =dk_local2_no_weekend, workplace ~ 
                  lag_workplace+ 
                  vote_share_populists+
                  after_march_10:vote_share_populists+
                  after_april_13:vote_share_populists+
                  after_may_10:vote_share_populists+
                  after_june_7:vote_share_populists+
                  log_house_prices+disp_inc+log_density+emp_finance_lib+emp_construction+emp_wholesale_retail_transport_acc+emp_adm_health_education+plus65+under18,
                fixed_effects = ~date+location, se_type = "HC1")

m3 <- lm_robust(data =dk_local2_no_weekend, workplace ~ 
                  lag_workplace+ 
                  after_march_10+after_april_13+after_may_10+after_june_7+
                  after_march_10:vote_share_populists+
                  after_april_13:vote_share_populists+
                  after_may_10:vote_share_populists+
                  after_june_7:vote_share_populists,
                fixed_effects = ~sub_region_2, se_type = "HC1")

m4 <- lm_robust(data =dk_local2_no_weekend, workplace ~ 
                  lag_workplace+ 
                  vote_share_populists:after_march_10+
                  after_april_13:vote_share_populists+
                  after_may_10:vote_share_populists+
                  after_june_7:vote_share_populists,
                fixed_effects = ~date+sub_region_2, se_type = "HC1")

m6 <- lm_robust(data =dk_local2_no_weekend, workplace ~ 
                  lag_workplace+
                  (vote_share_populists+log_house_prices+disp_inc+log_density+emp_finance_lib+emp_construction+emp_wholesale_retail_transport_acc+emp_adm_health_education+plus65+under18):after_march_10+
                  after_april_13:(vote_share_populists+log_house_prices+disp_inc+log_density+emp_finance_lib+emp_construction+emp_wholesale_retail_transport_acc+emp_adm_health_education+plus65+under18)+
                  after_may_10:(vote_share_populists+log_house_prices+disp_inc+log_density+emp_finance_lib+emp_construction+emp_wholesale_retail_transport_acc+emp_adm_health_education+plus65+under18)+
                  after_june_7:(vote_share_populists+log_house_prices+disp_inc+log_density+emp_finance_lib+emp_construction+emp_wholesale_retail_transport_acc+emp_adm_health_education+plus65+under18),
                fixed_effects = ~date+sub_region_2, se_type = "HC1")


coef_list <- list( "lag_workplace" = "Lag workplace" , 
                   "vote_share_populists" = "Populist vote" ,  
                   "disp_inc" = "Disposable income",
                   "log_house_prices" =  "Log (House prices)",
                   "emp_finance_lib" = "Employment: Financial services, real estate, liberal professions", 
                   "emp_construction" = "Employment: Construction",
                   "emp_wholesale_retail_transport_acc" = "Employment: wholesale, retail, transport, hotels, restaurants",
                   "emp_adm_health_education" = "Employment: defense, police, fire, education, health",
                   "under18" = "Pop < 18" ,  
                   "plus65" = "Pop > 65" , 
                   "log_density" = "Log (Density)" ,  
                   "after_march_10" = "After March 10",
                   "after_april_13" = "After April 13" , 
                   "after_may_10" = "After May 10" ,  
                   "after_june_7" = "After June 7",
                   "vote_share_populists:after_march_10" = "Populist vote * After March 10", 
                   "vote_share_populists:after_april_13" = "Populist vote * After April 14",
                   "vote_share_populists:after_may_10"  = "Populist vote * After May 10", 
                   "vote_share_populists:after_june_7" = "Populist vote * After June 7",
                   "after_march_10:vote_share_populists" = "Populist vote * After March 10", 
                   "after_april_13:vote_share_populists" = "Populist vote * After April 14",
                   "after_may_10:vote_share_populists"  = "Populist vote * After May 10", 
                   "after_june_7:vote_share_populists" = "Populist vote * After June 7")

# This table has Model 1 
#library(texreg)
texreg(list(m1, m2, m3, m4, m6), include.ci = FALSE, custom.coef.map = coef_list,
       custom.gof.rows = list("Region Effects"= c("Y", "Y", "N", "N", "N"), "Date Effects" = c( "N", "Y", "N","Y", "Y"),
                              "Municipality Effects" = c("N", "N", "Y", "Y", "Y"), "All Interactions" = c("N", "N", "N", "N", "Y"))
       , file="/Users/madselkjaer/Dropbox/WEALTHPOL_Research/Papers/WEP/figures/EU_analysis/DK_workplace.tex", caption = "Workplace Activity in Denmark", 
       caption.above = TRUE, label = "table_DK", booktabs = TRUE
)









#### APPENDIX MATERIAL: Adding positive test ratio to the models. 
# removing all NA's
dk_local2_no_na <- filter(dk_local2, pos_test_ratio!="NA") 

weekend_df_tests <- data.frame(date = unique(dk_local2$date)) %>%
  mutate(day = weekdays(date)) %>%
  mutate(weekend= case_when(
    day %in% c("Saturday", "Sunday") ~ 1,
    date==as.Date("2020-02-15") ~ 1,
    date==as.Date("2020-02-16") ~ 1,
    date==as.Date("2020-02-17") ~ 1,
    date==as.Date("2020-02-18") ~ 1,
    date==as.Date("2020-02-19") ~ 1,
    date==as.Date("2020-02-20") ~ 1,
    date==as.Date("2020-02-21") ~ 1,
    date==as.Date("2020-02-22") ~ 1,
    date==as.Date("2020-02-23") ~ 1,
    date==as.Date("2020-02-24") ~ 1,
    date==as.Date("2020-02-25") ~ 1,
    date==as.Date("2020-02-26") ~ 1,
    date==as.Date("2020-04-09") ~ 1,
    date==as.Date("2020-04-10") ~ 1,
    date==as.Date("2020-04-13") ~ 1,
    date==as.Date("2020-05-01") ~ 1,
    date==as.Date("2020-05-08") ~ 1,
    date==as.Date("2020-05-21") ~ 1,
    date==as.Date("2020-05-22") ~ 1,
    date==as.Date("2020-06-01") ~ 1,
    date==as.Date("2020-06-05") ~ 1,
    TRUE ~ 0))


## Figure A5
## coef plot by day
# workplace - DPP+NB+SK

coefs_map_workplace <- map_df(unique(dk_local2_no_na$date), function(date_i) {
  filtered_data <- filter(dk_local2_no_na, date == date_i) %>%
    lm(data =., workplace ~ vote_share_populists+log_house_prices+disp_inc+log_density+emp_finance_lib+emp_construction+emp_wholesale_retail_transport_acc+emp_adm_health_education+plus65+under18+pos_test_ratio+I(location)) %>%
    broom::tidy() %>%
    rename(coef = estimate, se = std.error) %>%
    mutate(date = date_i)
})

coefs_map_workplace <- coefs_map_workplace %>% 
  mutate(upper = coef+2*se,
         lower = coef-2*se) %>% 
  filter(term %in% c("vote_share_populists","log_house_prices", "disp_inc", "log_density", "under18", "pos_test_ratio")) %>% 
  mutate(term = case_when(term=="vote_share_populists" ~ "Populist Vote",
                          term=="log_house_prices" ~ "Log House Prices",
                          term=="disp_inc" ~ "Disposable Income",
                          term=="log_density" ~ "Log Density",
                          term == "under18" ~ "Population < 18",
                          term == "pos_test_ratio" ~ "Positive test ratio"))

coefs_map_workplace %>% 
  left_join(weekend_df_tests) %>%
  filter(weekend==0)  %>%
  ggplot(aes(x = date, y = coef))+
  geom_line()+
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
  geom_vline(xintercept =as.Date("2020-03-11"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-04-14"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-05-11"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-06-08"), linetype="dashed", color = "blue")+
  geom_ribbon(aes(ymin=lower, ymax = upper, alpha=0.3), color="lightgray")+
  #lockdown
  xlab("Date")+ylab("Day by Day Coefficients (DV: Workplace Activity)")+
  theme_minimal()+
  theme(legend.position = "none") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_wrap(~term, scales = "free_y")

